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Abstract 

We  study  a  population  model  in  which  there  are  two  species,  one  of  which  has  a 
juvenile  and  adult  life  stage.  The  adults  of  the  first  species  prey  on  the  second  species, 
which  in  turn  preys  on  the  juveniles  of  the  first.  One  version  of  the  model  represents 
systems  where  neither  species  can  survive  on  its  own,  although  we  find  that  both  can 
survive  through  mutual  predation.  To  avoid  extinction,  the  two  types  of  predation  must 
be  of  sufficient  strength  and  in  appropriate  proportion  to  one  another.  Another  version 
of  the  model  represents  systems  where  each  species  can  survive  without  the  other,  and 
there  we  find  that  mutual  predation  is  capable  of  increasing  both  of  their  equilibrium 
populations  or  creating  stable  limit  cycles. 

1  Introduction 

An  organism’s  trophic  level  is  the  position  that  it  occupies  in  the  food  web,  defined 
by  the  organisms  that  it  eats  and  vice  versa.  However,  few  organisms  occupy  the  same 
trophic  levels  throughout  their  lives.  Typically  an  organism’s  predators  and  prey  change 
over  its  lifetime,  as  do  the  organisms  of  similar  trophic  level  with  which  it  competes  for 
food.  Werner  and  Gilliam  [2]  review  many  examples  of  ecosystems  in  which  competitive 
and  predatory  relationships  are  age-  and  size-dependent.  As  an  instance  of  the  latter, 
adult  salamanders  and  newts  prey  on  one  another’s  juveniles.  Similarly,  frogs  eat  insects, 
while  insect  larva  eat  tadpoles.  Although  competition  is  important,  we  do  not  consider  it 
in  this  report,  where  our  object  is  to  explore  solely  the  dynamics  that  arise  in  ecological 
models  with  age-structured  predation. 


2  Development  of  the  models 

Any  model  for  age-structured  predation  must  have  at  least  two  species,  and  at  least 
one  of  the  species  must  have  multiple  life  stages.  Some  ecological  models  have  used 
continuous  age  or  size  variables  to  describe  life  stage,  which  requires  partial  differential 
equations  to  describe  the  change  of  populations  in  time,  but  the  simplest  possibility  is 
to  have  two  discrete  life  stages.  Thus,  the  simplest  possible  model  has  three  populations 
in  total:  one  species  with  one  life  stage  and  one  species  with  both  juvenile  and  adult 
life  stages.  For  concreteness  we  shall  call  our  populations  tadpoles  (T),  frogs  (F),  and 
insects  (/),  where  the  frogs  eat  the  insects,  the  insects  eat  the  tadpoles,  and  the  tadpoles 
and  frogs  beget  one  another  through  recruitment  and  reproduction.  For  simplicity  we 
shall  assume  that  recruitment  occurs  at  a  rate  proportional  to  tadpole  biomass,  that 
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reproduction  occurs  at  a  rate  proportional  to  frog  biomass,  and  that  a  constant  fraction 
of  the  biomass  being  transferred  through  predation  is  lost  to  metabolic  inefficiencies. 
The  model  as  described  thus  far  is  diagramed  in  Figure  1.  Through  the  “feeding/death” 
fluxes  in  Figure  1,  each  population  may  gain  biomass  by  feeding  on  the  external  envi¬ 
ronment  (which  excludes  the  other  two  populations  explicitly  modeled) ,  or  lose  biomass 
through  death.  All  that  remains  to  fully  define  our  population  model  is  to  specify  the 
functional  forms  of  the  predation  and  feeding/death  terms. 


Figure  1:  Biomass  flux  in  the  tadpole-frog-insect  system. 


2.1  Feeding  or  death  terms 

Assuming  the  tadpoles,  frogs,  and  insects  do  not  compete  with  one  another,  each  pop¬ 
ulation’s  feeding/death  term  should  be  independent  of  the  other  two  populations.  The 
simplest  choice  is  a  proportional  law,  which  in  the  absence  of  predation  creates  expo¬ 
nential  growth  or  decay,  depending  on  the  parameter  values.  Unbounded  exponential 
growth  is  unrealistic,  but  exponential  decay  is  feasible:  it  represents  a  case  where  the 
two  species  cannot  survive  without  predation,  and  we  analyze  such  a  model  in  Section  3. 
We  also  wish  to  consider  a  scenario  in  which  both  species  would  survive  without  the 
predation.  This  requires  that  the  system  have  carrying  capacities,  which  cause  initially 
exponential  growth  to  saturate  at  finite  values.  We  must  put  a  carrying  capacity  on  the 
insects,  but  we  have  a  choice  between  putting  one  on  the  frogs,  the  tadpoles,  or  both. 
We  shall  somewhat  arbitrarily  place  a  carrying  capacity  on  frogs  and  not  on  tadpoles. 
In  the  biological  reality  of  frogs  and  tadpoles  this  is  usually  accurate;  frogs  will  typically 
run  out  of  environmental  resources  before  tadpoles  do.  But  for  other  species,  resources 
may  certainly  be  scarcer  for  juveniles  than  for  adults.  We  analyze  the  carrying  capacity 
model  in  Section  4. 
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2.2  Predation  terms 


The  simplest  possible  predation  law  is  a  quadratic  one  in  which  the  rate  of  predation 
is  proportional  to  both  predator  and  prey  populations,  the  so-called  Holling  type  I 
functional  form  [1],  Some  models  of  age-structured  predation  employing  these  forms 
were  studied  by  Whitehead  and  Doering  [3],  but  such  a  form  fails  to  reflect  the  fact 
that  increasing  the  amount  of  prey  beyond  a  certain  point  stops  benefiting  the  predators 
because  there  is  a  limit  to  how  fast  they  can  eat  and  metabolize  the  prey.  The  Holling 
type  II  functional  form  avoids  this  problem  by  saturating  as  the  amount  of  prey  becomes 
large: 

„  ,  .  (prey  population)  (predator  population) 

rate  of  predation  =  C-i - - - - - —— - - — ; — - , 

1  +  C 2 (prey  population) 

where  C\  and  C2  are  constant  parameter.  We  have  examined  some  models  employing  the 
type  I  form,  but  it  leads  to  unbounded  growth  in  certain  cases,  so  we  shall  henceforth 
always  use  the  type  II  form.  It  is  not  hard  to  justify  needing  a  predation  law  with 
saturation,  but  it  is  certainly  not  clear  a  priori  that  the  type  II  form  is  the  best  choice. 
We  justify  this  choice  in  the  appendix,  where  we  derive  the  type  II  form  as  an  asymptotic 
limit  of  a  higher-dimensional  dynamical  system  and  numerically  compare  the  reduced 
system  with  the  higher-dimensional  system. 


3  Model  without  carrying  capacities 

Let  the  feeding/death  term  for  tadpoles  add  biomass  at  the  rate  exT,  where  ex  can  take 
either  sign,  and  likewise  for  frogs  and  insects.  Then,  the  dimensional  ODE  governing 
the  biomass  of  each  population  is 


IT 

T=-7T  +  /3F-k— - — 

1  +  l-iT 

(1) 

F  =  aT  (F  +  pp  A 

1  +  vl 

(2) 

■  IT  FI 

/  =  e,/  +  „«1  +  (i  T  \l  +  yV 

(3) 

where 


7  =  a  —  ex 
C  =  P  -  £f- 

Without  frogs  or  insects,  tadpoles  should  die  out,  and  likewise  for  frogs,  so  7  and  (  must 
be  positive.  Only  the  ej  parameter  may  be  negative.  If  a  <  7,  tadpole  feeding  is  adding 
biomass  to  the  system,  and  the  same  is  true  for  frogs  if  (3  <  (,.  The  parameters  ?y/  and  r/p 
are  metabolic  efficiencies,  so  they  must  fall  between  0  and  1.  We  now  nondimensionalize 
the  system  by 

t^H  TmIt  F 1— >  -3-F  (4) 

This  particular  nondimensionalization  is  attractive  because  it  retains  unique  coefficients 
on  all  the  nonlinear  predation  terms  and  minimizes  the  number  of  free  parameters 
(seven).  Its  drawback  is  that  the  different  populations  are  no  longer  in  equivalent  units 
of  biomass,  so  they  cannot  be  meaningfully  compared.  This  is  tolerable  because  we  are 
more  interested  in  qualitative  system  behavior  than  relative  population  biomasses.  The 
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nondimensional  ODE  is 
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3.1  Model  without  predation 

Without  predation,  the  model  is  linear,  and  the  T-F  system  decouples  from  the  insects: 


dL 

dt 


1  0\  /at\ 

-6  0  AF 
0  c/  \A// 


The  tadpoles  and  frogs  decay  exponentially  if  a  <  b  and  grow  exponentially  if  a  >  b, 
while  the  insects  decay  exponentially  if  c  <  0  and  grow  exponentially  if  c  >  0.  The 
behaviors  of  the  uncoupled  systems  split  the  full  system  rather  naturally  into  four 
separate  cases.  (Structurally  unstable  parameter  values  such  as  a  =  b  or  c  =  0  will 
not  be  considered  in  anything  that  follows.)  The  only  equilibrium  of  this  system  is  the 
origin,  and  the  decoupled  systems  are  rather  boring  without  predation.  Predation  can 
create  stable  fixed  points  and  limit  cycles  for  all  four  combinations  of  signs  of  (a  —  b) 
and  c,  but  only  the  case  in  which  the  T  —  F  and  I  systems  both  decay  without  predation 
is  biologically  reasonable,  so  we  shall  restrict  ourselves  to  this  case.  Henceforth  in  this 
section,  a  <  b  and  c  <  0. 

When  /  and  the  T-F  system  both  decay  without  predation,  the  mechanism  by 
which  predation  can  stabilize  the  system  is  the  following.  Suppose  that  biomass  enters 
the  system  though  frog  feeding  ((  <  (3)  and  leaves  the  system  through  tadpole  death 
(7  >  a),  with  tadpole  death  dominating  when  the  T-F  system  is  isolated.  In  other 
words,  the  rate  of  reproduction  is  higher  than  ideal.  Introducing  insects  creates  a  flux 
of  biomass  from  tadpoles  to  frogs,  where  it  is  used  to  create  more  biomass  from  external 
feeding,  rather  than  lost  to  tadpole  death.  This  effect  can  prevent  the  entire  system 
from  decaying  to  zero,  even  with  the  additional  sources  of  biomass  loss  by  insect  death 
and  metabolic  inefficiencies. 


3.2  Lyapunov  bound 

Not  only  is  the  origin  linearly  stable  when  a  <  b  and  c  <  0,  the  Lyapunov  functional 
L  =  |F+^F  +  /  suffices  to  show  that  all  solutions  decay  to  the  origin.  This  is  a  rather 
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narrow  range  of  validity,  but  it  suggests  the  possible  importance  of  the  ratio  In 
terms  of  dimensional  variables, 

eg  f3 


"77  =  VIV  F~- 
df  7 

Biologically,  gigF  is  the  squared  geometric  mean  of  the  metabolic  efficiencies,  so  we  can 
think  of  this  term  roughly  as  the  metabolic  efficiency  of  the  entire  system. 


3.3  Nontrivial  equilibria 


At  nontrivial  fixed  points, 


F  =  T 


dl 

1  +  T 


aT  =  F 


0  =  c  + 


eT 

1  +  T 


fF 
1  +  /' 


(8) 

(9) 

(10) 


From  Equation  (8),  I  and  T  uniquely  define  F  at  a  fixed  point.  Using  this  equation  to 
eliminate  F  from  the  latter  two  yields  two  polynomials  equations  in  /  and  T: 


o(l  +  I)(l  +  T)  =  (1  +  T  +  dl){b  +bl-  gl)  (11) 

c(l  +  I)  (1  +  T)  +  eT(l  +  J)  =  /T(l  +  T  +  dl) .  (12) 

Equation  (11)  is  linear  in  T,  and  Equation  (12)  is  linear  in  /,  so  they  may  be  used  to 
find  explicit  expressions  for  T(T)  and  /(T),  respectively: 

T_  dl(b+(b-g)l)  1 
(a  +  g  —  b)I  +  a  —  b 
7_  fT2  —  (c  +  e  —  f)T  —  c . 

(c  +  e  —  df)T  +  c 

Evidently,  /  and  T  define  one  another  uniquely.  Applying  T{I)  to  Equation  12  gives 
a  cubic  polynomial  for  /  at  the  equilibria,  so  there  are  at  most  three  positive  real 
population  equilibira. 

i3I 3  +  i2I2  +  ill  +  io  =  0,  where 


10  =  —  e(a  —  b)2  <  0 

11  =  (a  —  b)  [2  e(b  —  g)  +  bd(c  +  e)  +  a(df  —  2e)] 

%2  =  d(c  +  e)  [(a  —  b)(b  —  g)  +  b{a  +  g  —  b)]  —  e(a  +  g  —  b)2  +  adf{a  +  g  —  b  —  bd ) 
i3  =  d(b  -  g)[(c  +  e)(a  +  g  -  b)  -  adf } 

The  cubic  equation  for  I  of  course  has  explicit  algebraic  solutions,  but  they  are  too 
messy  to  be  of  use,  so  although  we  have  explicit  expressions  for  F  and  T  in  terms  of 
/,  we  cannot  generally  predict  whether  F  and  T  will  be  positive  when  I  is  positive. 
The  best  we  can  do  analytically  is  infer  some  partial  information  about  the  signs  of 
polynomial  roots.  For  this,  we  also  need  the  first  and  last  coefficients  of  the  cubic 
polynomial  governing  T,  which  we  derive  by  applying  7(T)  to  Equation  11: 

t3T3  +  t2T2  +  t\T  +  T  =  0,  where 


T  =  —c2g  <  0 

h  =  f[(c+  e)(a  +  g  -  b)  -  adf}. 
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Theorem  Let  a  <b  and  c  <  0.  If  g  <  b  and  adf  <  (c  +  e)(a  +  g  —  b),  there  are  either 
1  or  3  positive  equilibria.  Otherwise,  there  are  either  0  or  2  positive  equilibria. 

Proof  The  I  polynomial  is  cubic,  and  I  uniquely  defines  F  and  T,  so  there  are  at 
most  three  real,  positive  equilibria.  It  is  clear  from  Equation  8  that  F  is  positive  when 
I  and  T  are,  so  it  suffices  to  know  when  the  I  and  T  polynomials  have  corresponding 
positive  roots.  There  are  three  ways  in  which  the  number  of  positive  roots  of  the  I  or  T 
polynomials  may  change.  Firstly,  a  root  may  remain  real  but  leave  the  positive  octant 
if  ?'o  or  to  becomes  negative,  but  we  have  eliminated  this  possibility  by  assumption. 
Secondly,  the  number  of  positive  equilibria  may  change  by  two  when  a  pair  of  equilibria 
become  complex  simultaneously,  a  saddle-node  bifurcation.  Thirdly,  an  equilibrium 
may  move  off  to  infinity  when  i3  or  f3  pass  through  zero.  By  dividing  parameter  space 
into  four  regions  in  which  i3  and  to  do  not  change  sign,  we  are  assured  that  the  number 
of  equilibria  within  each  region  may  change  only  by  saddle-node  bifurcations.  These 
regions  are 


I  =  {i 3  >  0,  to  >  0} 
II  =  {?3  <0,  f3  >  0} 

III  =  {*3  >0,  to  <  0} 

IV  =  {*3  <  0,  t3  <  0} 


{g  <  6,  adf  <  (c  +  e)(a  +  g  -  b)} 

{g  >  6,  adf  <  (c  +  e)(a  +  g  -  b)} 

{g  <  6,  adf  >  (c+  e)(a  +  g  -  b)} 

{g  >  6,  adf  >  (c+  e)(a  +  g  -  b)}. 


If  we  know  the  number  of  positive  equilibria  at  one  point  in  each  parameter  region,  we 
know  the  number  of  positive  equilibria  at  all  points  in  that  parameter  region,  modulo 
2.  So,  we  chose  one  such  test  point  in  each  parameter  region  and  computed  the  equi¬ 
libria  numerically,  inferring  from  this  the  possible  number  of  positive  equilibria  in  each 
parameter  region.  The  results  are  tabulated  in  Table  1.  Evidently,  there  may  be  1  or 
3  positive  equilibria  in  parameter  region  I,  and  0  or  2  otherwise.  Parameter  region  I 
is  defined  precisely  by  the  condition  that  g  <  b  and  adf  <  (c  +  e)  (a  +  g  —  b) ,  so  the 
theorem  is  proved.  I 


Table  1:  Number  of  nontrivial  equilibria  in  each  of  the  four  parameter  regions,  as  inferred 
from  a  test  points,  (a,  b,  c,  d,  e,  /,  g),  in  parameter  space. 


Parameter  region 

Test  point 

Number  of  positive 
equilibria  at  test  point 

Number  of  positive 
equilibria  in  region 

I 

(0.9,1, -0.5, 1,1.5, 0.1, 0.5) 

1 

1  or  3 

II 

(0.9,1, -0.5, 1,1.5, 0.1, 2) 

2 

0  or  2 

III 

(0.9,1,-0.5,100,1.5,0.1,0.5) 

0 

0  or  2 

IV 

(0.9,1,-0.5,100,1.5,0.1,2) 

0 

0  or  2 

By  the  preceding  theorem,  we  have  partial  information  about  the  number  of  equilib¬ 
ria  inside  and  outside  region  I.  Because  a  saddle-node  bifurcation  creates  two  equilibria 
of  opposite  stability,  we  also  know  that  there  can  be  at  most  two  stable  equilibria  inside 
parameter  region  I,  and  at  most  one  stable  equilibrium  outside  it.  We  would  of  course 
like  precise  conditions  on  when  the  saddle-node  bifurcations  occur,  thereby  further  di¬ 
viding  our  parameter  regions  into  ones  in  which  the  numbers  of  positive  equilibria  are 
exactly  known.  This  requires  knowing  when  I  or  T  become  complex.  Since  I  and  T  are 
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governed  by  cubic  equations,  this  can  in  principle  be  determined  from  the  discriminants 
of  those  equations.  The  I  equation  has  one  real  and  two  complex  roots  precisely  when 

A/  =  18*3*2*i*o  —  4»2?o  +  i\i\  —  4 *3*f  —  27*3*q  <  0, 

and  likewise  for  the  coefficients  of  the  T  equation,  so  the  parameter  regimes  we  seek 
are  divided  by  the  surfaces  on  which  the  discriminants  vanish.  Unfortunately,  these 
inequalities  are  prohibitively  messy  when  expressed  in  terms  of  the  problem  parameters. 
In  biological  modeling,  one  is  more  interested  in  qualitative  behaviors  than  the  precise 
values  at  which  bifurcations  occur,  so  we  will  be  content  to  observe  the  saddle-node 
bifurcations  numerically. 

3.4  Numerical  exploration  of  parameter  space 

Asymptotic  analysis  shows  that  the  predation  parameters  cannot  be  small  compared  to 
the  population  decay  rates  if  predation  is  to  stabilize  the  decay.  In  fact,  numerical  exper¬ 
iments  reveal  that  the  predation  parameters  must  typically  be  about  an  order  of  mag¬ 
nitude  larger  than  the  decay  rates.  In  light  of  this,  we  shall  fix  (a,b,c)  =  (0.9,1,— 0.1) 
and  expect  interesting  behavior  when  predation  parameters  are  0(1). 

3.4.1  Looking  for  bifurcations 

When  ( d,e,f,g )  =  (0.2, 1,  0.2, 0.5),  there  are  three  positive  equilibria,  the  maximum 
number  possible,  so  we  choose  this  as  a  starting  point  from  which  to  explore  the  four¬ 
dimensional  space  of  predation  parameters.  Linearizing  the  system  about  these  equilib¬ 
ria  and  computing  matrix  eigenvalues  numerically,  we  find  that  none  of  the  fixed  points 
are  stable.  Biologically,  we  are  interested  in  regimes  where  there  are  stable  fixed  points, 
so  we  use  the  bifurcation  continuation  package  MATCONT  to  continue  these  equilibria 
in  parameter  space.  Arbitrarily  choosing  the  parameter  e  in  which  to  continue  the  equi¬ 
libria,  we  obtain  the  bifurcation  diagram  of  Figure  2.  We  have  chosen  the  coordinate  I 
for  the  ordinate  of  our  bifurcation  diagrams  because  the  value  of  /  sometimes  becomes 
unrealistically  small,  and  we  wish  to  see  when  this  is  so. 

Continuing  in  e,  we  find  a  saddle-node  bifurcation  ( LP ),  a  subcritical  Hopf  bifurca¬ 
tion  ( H+ ),  and  a  neutral  saddle  (NS).  These  are  all  bifurcations  of  codimension  1,  as 
we  would  expect  to  find  when  varying  only  one  parameter.  To  access  the  higher-level 
structure  of  parameter  space,  we  would  like  to  find  higher-codimension  bifurcations. 
The  MATCONT  package  can  in  general  only  find  bifurcations  up  to  codimension-2,  so 
we  shall  settle  for  this,  though  ideally  we  would  like  to  find  bifurcations  of  codimension 
up  to  the  dimension  of  our  parameter  space.  Codimension-2  bifurcations  are  typically 
found  by  continuing  a  codimension- 1  bifurcation  in  two  parameters,  so  we  shall  continue 
all  of  the  bifurcations  of  Figure  2.  Continuing  in  e  and  g  yields  the  bifurcation  diagram 
of  Figure  3,  in  which  we  see  two  types  of  codimension-2  bifurcations:  Bogdanov-Takens 
bifurcations  ( BT ),  and  generalized  Hopf  bifurcations  (GH),  also  known  as  Bautin  bifur¬ 
cations.  Continuing  the  codimension- 1  bifurcations  of  Figure  2  in  any  other  combination 
of  two  parameters  does  not  yield  any  other  types  of  codimension-2  bifurcations. 

In  the  neighborhood  of  a  Bogdanov-Takens  bifurcation,  there  are  guaranteed  to  be 
a  saddle-node  bifurcation  and  a  Hopf  bifurcation  (both  of  which  we  have  seen  already), 
and  also  a  saddle  homoclinic  bifurcation.  In  the  neighborhood  of  a  generalized  Hopf 
bifurcation,  there  are  guaranteed  to  be  both  supercritical-  and  subcritical  Hopf  bifur¬ 
cations,  and  a  fold  bifurcation  of  limit  cycles.  Although  we  are  only  assured  of  these 
system  behaviors  in  local  neighborhoods  of  the  codimension-2  bifurcations,  we  can  rea¬ 
sonably  expect  to  see  them  all  around  parameter  space. 
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Figure  2:  Continuation  in  e  of  the  three  starting  equilibria. 


Figure  3:  Continuation  in  e  and  g  of  the  Hopf  bifurcation  of  Figure  2. 
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Before  exploring  the  effects  of  each  parameter,  we  would  like  to  find  a  baseline  point 
in  parameter  space  at  which  all  parameters  are  of  similar  magnitude,  and  at  which  a 
stable  fixed  point  exists.  There  is  always  a  stable  fixed  point  near  a  Hopf  bifurcation, 
so  we  continue  the  Hopf  bifurcation  of  Figure  2  in  various  pairings  of  the  predation 
parameters,  arriving  at  a  point  where  (d,e,f,g)  is  near  (2, 4, 2,  2).  We  shall  use  this 
and  other  nearby  values  as  our  baseline  parameter  points. 

3.4.2  Overall  predation  strength 

To  study  the  effect  of  the  overall  predation  strength  in  the  system,  we  fix  the  ratios 
between  all  four  predation  parameters  and  vary  them  proportionally.  The  resulting 
diagrams  are  shown  in  Figure  4,  where  we  have  used  two  sightly  different  parameter 
ratios  to  exhibit  one  case  each  where  the  Hopf  bifurcation  is  supercritical  or  subcritical. 
In  any  of  the  bifurcation  diagrams  that  follow,  the  predation  parameter  ratios  could 
be  chosen  to  realize  either  type  of  Hopf  bifurcation,  but  the  distinction  is  not  very 
important  because  limit  cycles  in  this  system  are  not  robust  biologically,  for  the  following 
reason.  In  the  chosen  regime  of  a,  b  and  c,  the  coordinates  of  the  stable  fixed  point 
appear  to  always  be  such  that  I  <C  T,F.  As  the  limit  cycle  emerging  from  a  Hopf 
bifurcation  grows,  the  oscillations  become  strongly  nonlinear,  consisting  of  a  fast  part 
and  slow  part.  During  the  slow  part  of  the  cycle,  the  value  of  I  is  very  small.  Without 
traveling  very  far  in  parameter  space,  either  the  limit  cycle  is  destroyed  in  a  homoclinic 
bifurcation,  or  /  becomes  so  small  during  the  slow  part  of  the  cycle  that  in  the  non¬ 
continuum  of  reality,  the  insect  population  would  go  extinct,  bringing  the  same  fate 
to  the  frogs  and  tadpoles.  If  we  altered  the  model  to  include  stochasticity  or  discrete 
populations,  we  could  reasonably  expect  it  to  display  this  biologically-realistic  extinction 
during  such  troughs  of  I.  The  type  of  Hopf  bifurcation  controls  whether  or  not  a  stable 
limit  cycle  coexists  with  a  stable  fixed  point,  but  since  the  limit  cycle  is  a  much  less 
robust  structure  than  the  fixed  point,  the  distinction  is  not  very  important. 

It  is  evident  from  Figure  4  that  the  system  has  no  stable  fixed  points  or  limit 
cycles  when  the  predation  is  too  weak.  This  is  expected,  since  the  predation  must 
stabilize  the  exponential  decays  of  the  uncoupled  systems,  and  it  makes  clear  that 
predation  cannot  do  this  as  a  mere  perturbation  on  the  uncoupled  systems.  As  predation 
increases,  the  fixed  point  not  only  remains  stable,  but  its  basin  of  attraction  enlarges. 
Simultaneously,  however,  the  equilibrium  insect  population  decreases  toward  zero,  and 
the  frog  and  tadpole  populations  decrease  asymptotically  toward  identical  finite  values. 
When  the  basin  of  attraction  is  too  small,  a  perturbation  could  send  the  system  across 
the  separatrix  and  onto  a  trajectory  heading  toward  the  origin  (extinction).  When 
the  equilibrium  insect  population  is  too  small,  it  is  vulnerable  to  eradication  by  some 
catastrophic  event,  in  which  fate  the  frogs  and  tadpoles  would  follow.  The  fitness  of  the 
insect  population  would  be  maximized  at  some  intermediate  predation  strength  that 
balances  these  two  factors. 

3.4.3  Relative  strengths  of  the  two  types  of  predation 

The  parameters  d  and  e  convey  the  strength  of  the  insect-tadpole  predation,  and  the 
parameters  /  and  g  do  likewise  for  the  frog-insect  predation.  To  study  the  effect  of  the 
relative  strengths  of  the  two  types  of  predation,  we  vary  the  insect-tadpole  predation 
parameters  in  fixed  proportion  while  holding  the  other  two  constant,  and  likewise  for  the 
frog-insect  parameters.  The  resulting  bifurcation  diagrams  are  presented  in  Figure  5, 
and  they  both  suggest  the  same  conclusions.  The  process  in  which  insects  move  biomass 
from  tadpoles  to  frogs  is  essential  to  stabilizing  the  system,  and  when  I-T  predation 
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Figure  4:  Bifurcation  diagrams  in  which  all  four  predation  parameters  are  varied  propor¬ 
tionally  to  study  the  overall  effect  or  predation  strength.  Different  proportions  between  the 
parameters  can  produce  either  supercritical  (top)  or  subcritical  (bottom)  Hopf  bifurcations. 
Light  gray  depicts  unstable  limit  cycles,  while  dark  gray  depicts  stable  limit  cycles. 
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is  too  weak  relative  to  F-I  predation,  the  insect  population  is  kept  low  and  cannot 
transfer  enough  biomass,  so  the  stable  fixed  point  disappears  through  a  saddle-node 
bifurcation.  Conversely,  when  the  I-T  predation  is  relatively  strong,  the  fixed  point 
will  lose  stability  to  stable  oscillations.  In  such  oscillations,  insects  eat  tadpoles  quickly 
and  deplete  the  tadpole  population,  after  which  frogs  eat  the  plentiful  insects,  depleting 
that  population,  and  then  frogs  give  birth  to  many  tadpoles,  losing  biomass  themselves 
until  the  insects  rebound.  When  these  oscillations  become  too  dramatic,  extinction  is 
the  likely  result. 

4  Model  with  carrying  capacities 

With  carrying  capacities  on  frogs  and  insects,  the  dimensional  governing  equations 
become 

T  =  -7  T  +  f3F  —  k 

1  +  /iT 

FI 

F=aT-(F(l+jrF)+VF\YT-i 

IT  FI 

1  =  'U(1  -  ip)  +  ""‘nyr  ~  XT+Vr 

4.1  Nontrivial  equilibria 

Without  predation,  a  nontrivial  equilibrium  exists  only  if  ^  >  1,  which  is  also  the 
condition  under  which  the  origin  is  unstable.  We  wish  to  study  a  parameter  regime  in 
which  the  system  has  nontrivial  behavior  before  predation  is  added,  so  we  shall  assume 
this  inequality  always  holds.  With  predation,  a  nontrivial  fixed  point  must  satisfy 

F=1Fb+^) 

These  algebraic  equations  cannot  be  solved  explicitly  for  the  equilibrium  populations, 
so  we  must  resort  to  perturbation  expansions  and  numerical  solutions. 

4.1.1  Perturbation  of  the  nontrivial  equilibrium  by  weak  predation 

Let  the  predation  strength  be  small: 

k  =  £K\  A  =  eAi, 

where  eCl.  We  can  expand  the  equilibrium  populations  in  e  (e.g.  T  ~  To+eTi+0(e2)), 
going  to  first  order  to  obtain  the  leading  order  impact  of  predation.  At  zeroth  order, 
i.e.  without  predation, 


T0=  ?(s£-l)JV 

(13) 

*b=($  ~l)N 

(14) 

h  =  M. 

(15) 
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Figure  5:  Bifurcation  diagrams  in  which  the  relative  strengths  of  the  two  types  of  predation 
are  varied.  To  study  the  variation  of  insect-tadpole  predation  strength  (top),  d  and  e  are 
varied  proportionally  (e  =  2d),  while  frog-insect  predation  is  kept  constant  (/  =  2,  g  =  2.5). 
To  study  the  variation  of  frog-insect  predation  (bottom),  /  and  g  are  varied  proportionally 
{f  =  \g),  while  insect-tadpole  predation  is  kept  constant  (d  =  2,  e  =  4).  Light  gray  depicts 
unstable  limit  cycles,  while  dark  gray  depicts  stable  limit  cycles. 
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The  first  order  coefficients  in  the  expansions  are 
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Unlike  the  system  without  carrying  capacities,  this  system  has  stable  fixed  points 
without  predation,  so  we  may  compare  populations  with  and  without  predation.  There 
are  parameter  regimes  in  which  only  amphibians,  or  only  insects,  benefit  from  the 
predation.  There  is  also  a  regime  in  which  predation  hurts  both  species,  for  instance 
when  both  species  are  bad  at  metabolizing  each  other.  This  situation  is  essentially 
a  prisoner’s  dilemma;  both  species  would  benefit  from  a  truce,  but  once  a  system  of 
reciprocal  predation  has  evolved,  it  is  hard  to  see  how  it  could  stop.  Finally,  there  is  a 
parameter  regime  where  the  predation  benefits  both  species  by  improving  the  efficiency 
with  which  the  entire  tadpole-frog-insect  system  uses  environmental  resources.  We  will 
call  this  a  regime  of  population  increase  by  mutual  predation,  or  PIMP,  and  we  now 
look  at  this  regime  in  further  detail. 


4.1.2  Population  increase  by  mutual  predation 

Suppose  the  two  types  of  predation  occur  with  relative  strengths  given  by  K  =  A .  We 
consider  turning  on  predation  with  fixed  K,  so  the  signs  of  the  total  derivatives  of  T, 
F ,  and  I  with  respect  to  k  (or  A)  dictate  whether  predation  increases  or  decreases  the 
population  biomasses  at  equilibrium.  When  e  is  small, 
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and  likewise  for  the  other  populations.  Evaluated  at  k  =  0,  just  as  predation  is  turned 
on, 
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The  signs  of  the  above  derivatives  determine  whether  predation  hurts  or  helps  a  species, 
insofar  as  its  equilibrium  biomass  decreases  or  increases.  The  regime  of  PIMP  could  be 
defined  by  the  condition  that  both  species’  biomasses  increase  with  predation,  or  by  the 
stronger  condition  that  all  three  populations’  biomasses  increase,  and  it  is  not  obvious 
which  definition  is  more  useful.  Biologically,  the  advantage  of  a  larger  population  is 
that  it  is  more  fit  because  it  is  more  genetically  diverse  and  robust  to  catastrophe.  In 
that  sense,  it  does  not  matter  what  life  stage  the  amphibians  are  in  because  one  life 
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stage  begets  another,  and  our  two-stage  model  is  a  simplification  anyway.  On  the  other 
hand,  it  could  still  be  catastrophic  for  an  entire  generation  of  tadpoles  to  be  wiped  out. 
We  will  thus  derive  the  conditions  for  PIMP  using  both  the  strong  and  weak  definitions. 


4.1.3  Conditions  for  strong  population  increase  by  mutual  predation 

All  three  derivatives  with  respect  to  n  are  positive,  meaning  the  equilibrium  biomasses 
of  insects  and  both  amphibian  life  stages  increase,  if  and  only  if 
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(16) 


The  ratio  of  predation  rates,  K ,  and  the  rate  of  reproduction,  (3,  can  be  changed  by  be¬ 
havior,  so  they  can  vary  on  much  shorter  time  scales  than  the  other  parameters,  which 
are  controlled  by  physiology.  Thus,  we  may  regard  K  and  (3  as  control  parameters.  Re¬ 
calling  that  (  =  0  —  6t,  it  is  clear  from  Equations  (13)  and  (14)  that  without  predation, 
changing  the  birth  rate  cannot  increase  frog  and  tadpole  populations  simultaneously; 
only  predation  can  do  that.  There  exists  an  interval  of  K  in  which  strong  PIMP  occurs 
precisely  when 

^(2^_1) <mF-  (17) 

For  this  inequality  to  hold,  it  is  necessary  that  metabolic  losses  not  be  too  large,  and 
also  that  a  <  7  and  0  >  That  is,  predation  can  only  increase  all  three  equilibrium 
biomasses  when  frogs  gain  biomass  through  their  interaction  with  the  external  environ¬ 
ment,  while  tadpoles  lose  it.  The  biological  interpretation  of  this  parameter  regime  is 
that  frogs  are  more  fit  for  their  environment  than  tadpoles.  Indeed,  one  could  expect 
this  to  be  true;  there  are  morphogenetic  tradeoffs  between  juvenile  and  adult  fitness, 
and  since  frogs  spend  the  majority  of  their  lives  in  the  adult  stage,  it  is  likely  that  they 
would  evolve  to  be  maximally  fit  as  adults.  It  is  less  clear  that  tadpoles  would  be  so 
unfit  as  to  lose  biomass  without  a  constant  input  from  reproduction,  but  this  is  certainly 
feasible  in  harsh  environments  that  create  high  juvenile  mortality.  So,  in  such  a  regime, 
predation  can  increase  all  three  biomasses  by  increasing  the  ratio  of  frog  biomass  to 
tadpole  biomass,  with  insects  profiting  as  middlemen. 


4.1.4  Conditions  for  weak  population  increase  by  mutual  predation 

Let  A  =  T  +  F,  the  total  amphibian  biomass.  Without  predation,  the  equilibrium 
amphibian  biomass  is 
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There  is  an  optimal  birth  rate,  0 ,  that  maximizes  Aq  .  In  the  strong  definition  of  PIMP 
there  is  no  such  way  to  define  an  optimal  0  since  there  is  always  a  tradeoff  between  T 
and  F.  When  /3  is  above  its  optimal  value,  too  much  of  the  frogs’  biomass  is  going  into 
tadpoles,  who  are  less  fit  than  the  frogs.  When  (3  is  below  its  optimal  value,  the  frog 
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population  becomes  so  large  than  environmental  pressures  make  frogs  less  successful 
than  tadpoles. 

The  dependence  of  equilibrium  amphibian  biomass  on  weak  predation  is 
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Weak  predation  increases  the  equilibrium  values  of  both  A  and  I,  thus  satisfying  the 
definition  of  weak  PIMP,  if  and  only  if 
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and  clearly  such  a  K  exists  if  and  only  if 
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This  condition  is  harder  to  interpret  biologically  than  the  stronger  condition  of  Equa¬ 
tion  (eq:  strong  PIMP  necessary) ,  but  it  has  roughly  the  same  necessary  conditions;  in 
the  parameter  regime  for  which  predation  improves  efficiency,  the  metabolisms  cannot 
be  too  inefficient,  frogs  cannot  be  too  unfit,  and  tadpoles  cannot  be  too  fit. 

Adding  predation  is  like  decreasing  the  rate  of  reproduction  in  that  it  transfers 
biomass  from  tadpoles  to  frogs,  though  it  is  certainly  not  identical.  We  have  not  fully 
explored  the  relationship  between  optimal  /3  and  optimal  A',  but  this  would  be  a  good 
topic  for  future  work.  For  instance,  we  would  like  to  know  whether  Equation  (18) 
can  hold  when  the  reproduction  rate  is  at  its  no-predation  optimal  value,  and  how  the 
optimal  reproduction  rate  changes  in  the  presence  of  predation. 


4.1.5  Beyond  weak  predation 


We  have  seen  analytically  that  there  is  a  parameter  regime  in  which  PIMP  occurs  when 
predation  is  weak,  but  we  would  like  to  know  whether  increasing  the  strength  of  the 
predation  will  increase  the  equilibrium  populations  indefinitely.  We  have  not  tackled  this 
question  analytically,  but  in  all  numerical  experiments  the  populations  reach  a  maximum 
value  before  decreasing  as  predation  is  strengthened  further.  To  use  MATCONT,  we 
prefer  to  work  with  dimensionless  parameters,  so  we  apply  the  nondimensionalization 
of  Equation  (4),  the  same  one  used  for  the  system  without  carrying  capacities. 
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TT  FT 
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where  M  and  N  have  been  nondimensionalized  in  the  same  way  as  I  and  F,  respectively. 

In  terms  of  dimensionless  quantities,  the  regime  of  strong  PIMP  given  by  Equa¬ 
tion  (16)  becomes 


(2a 


,,d 

V 


l  +  (f  -1)« 

1  +  M 


An  example  case  that  falls  in  this  regime  has  (a,  b,  c )  = 
predation  parameters  in  fixed  proportion  such  that  e  = 


e 


(2,1,1),  (M,  N)  = 
-  2d,  /  =  d,  and  g 


(1,1),  and  the 
=  4d.  Figure  6 


227 


displays  the  variation  of  equilibrium  populations  with  predation  strength.  Starting 
with  no  predation  and  increasing  the  predation  parameters  in  fixed  proportion,  we 
initially  see  increases  in  all  three  equilibrium  populations,  as  predicted  by  the  asymptotic 
analysis,  followed  by  decreases  in  all  three  populations.  Insect  population  appears  to 
go  to  zero  as  predation  goes  to  infinity,  while  frog  and  tadpole  populations  appear  to 
decrease  asymptotically  to  the  same  finite  value.  Although  it  might  be  possible  for  the 
amphibians  to  predate  the  insects  to  extinction,  it  is  not  in  their  interest  to  do  so;  this 
strategy  does  not  maximize  their  equilibrium  population. 


Figure  6:  Equilibrium  populations  as  predation  strength  is  varied  in  the  strong  PIMP 
regime.  Fixed  parameters  parameters  are  (a,  b,  c )  =  (2, 1, 1)  and  (M,  N)  =  (1, 1).  Predation 
parameters  are  varied  in  constant  proportion  with  e  =  2d,  f  =  d,  and  g  =  4 d. 


4.2  Numerical  exploration  of  parameter  space 

We  have  searched  for  bifurcations  in  parameter  space  using  MATCONT,  but  we  have  not 
found  any  in  regimes  where  PIMP  occurs;  we  simply  observe  the  equilibrium  changing 
its  coordinates  while  remaining  stable.  In  other  parameter  regimes  we  find  supercritical 
Hopf  bifurcations  but  no  saddle-node  bifurcations.  The  Hopf  bifurcations  observed  in 
this  system  were  of  a  different  character  than  those  in  the  system  without  carrying 
capacities.  In  that  system,  continuing  in  parameter  space  past  a  Hopf  bifurcation  typi¬ 
cally  led  to  growing  oscillations  and,  ultimately,  extinction.  In  the  system  with  carrying 
capacities,  a  limit  cycle  may  appear  and  grow  as  we  move  through  a  Hopf  bifurcation, 
but  as  we  continue  through  parameter  space,  it  ultimately  shrinks  and  disappears  back 
into  a  stable  equilibrium.  So,  it  seems  limit  cycles  occur  only  in  isolated  regions  of 
parameter  space  where  certain  resonances  are  strong.  This  significant  qualitative  dif- 
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ference  from  the  system  without  carrying  capacities  is  due  to  the  damping  effect  of  the 
carrying  capacities. 


5  Conclusions 

In  the  model  without  carrying  capacities,  we  examined  the  regime  in  which  the  insect 
and  frog-tadpole  systems  decay  without  predation.  Predation  can  create  stable  fixed 
points  and  limit  cycles,  though  the  limit  cycles  are  not  robust  and  rather  unrealistic 
biologically,  so  in  this  model  the  two  species  are  completely  codependent.  For  this 
stabilization  to  occur,  predation  must  be  above  a  certain  minimum  strength,  and  the 
relative  strengths  of  insect-tadpole  and  frog-insect  predation  must  be  within  a  certain 
interval.  There  are  biological  parameters  for  which  a  fixed  point  and  a  limit  cycle  are 
simultaneously  stable,  but  we  have  not  found  parameters  for  which  nontrivial  stable 
fixed  points  coexist. 

In  the  model  with  carrying  capacities,  both  species  exist  stably  without  predation. 
Depending  on  the  biological  parameters,  predation  may  increase  or  decrease  the  equi¬ 
librium  populations  in  any  combination.  Fixing  all  biological  parameters  other  than 
predation  rates,  we  have  seen  that  mutual  predation  can  increase  all  of  the  equilibrium 
populations  simultaneously.  However,  there  always  exists  an  optimal  rate  of  reproduc¬ 
tion  that  maximizes  the  equilibrium  biomass  of  amphibians,  and  we  have  not  determined 
whether  predation  can  only  increase  all  populations  simultaneously  by  effectively  de¬ 
creasing  the  rate  of  reproduction  toward  its  optimal  value.  If  this  is  the  case,  we  must 
ask  whether  it  is  biologically  realistic  for  the  amphibians  to  habitually  reproduce  at  a 
non-optimal  rate,  and  if  not,  we  would  conclude  that  population  increase  by  mutual 
predation  is  not  biologically  relevant  in  the  long  term.  This  is  both  a  mathematical  and 
biological  question  for  future  study. 
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7  Appendix:  The  Holling  type  II  functional  form 

The  simplest  system  in  which  the  Holling  type  II  functional  form  could  appear  is  a  2D 
predator-prey  model.  However,  the  justification  for  using  this  specific  form  in  such  a 
model  is  not  clear  a  priori.  A  quadratic  predation  term  that  is  proportional  to  both 
predator  and  prey  populations  would  not  require  much  justification  since  it  represents 
leading  order  behavior  at  the  very  least,  but  such  a  term  produces  unrealistic  exponential 
growth  in  the  system.  However,  building  on  the  work  of  Whitehead  and  Doering  [4],  we 
can  create  a  dichotomy  between  hungry  predators  sated  predators  and  use  only  linear 
and  quadratic  terms  in  the  resulting  3D  system,  and  then  find  that  in  the  appropriate 
limit  it  reduces  to  a  2D  predator-prey  system  with  Holling  type  II  predation  laws.  The 
functional  forms  of  the  3D  model  do  not  require  justification  since  they  are  the  simplest 
laws  that  could  capture  the  necessary  behavior.  To  evaluate  the  accuracy  of  the  reduced 
system  with  Holling  type  II  laws,  we  shall  compare  its  features  to  those  of  the  full  3D 
system. 
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7.1  The  3D  system 

The  3D  system  consists  of  hungry  predator  ( H ),  sated  predators  (S),  and  prey  (P): 

H  =  —(f>PH  +  gS  -  SH  +  (3S 
S  =  (j)PH  -gS-dS 
P=-0PJf  +  7P(  l-£). 

Hungry  predators  eat  prey  at  rate  (f>  and  become  sated  predators,  while  sated  predators 
metabolize  and  become  hungry  at  rate  fi.  Both  hungry  and  sated  predators  die  at  rate 
6,  and  sated  predators  give  birth  to  hungry  ones  at  rate  (3  >  5.  The  prey  is  being 
born  according  to  a  logistic  law.  In  a  more  accurate  model,  the  (j>  and  5  constants  in 
the  different  equations  could  take  different  values,  but  this  would  only  add  parameters 
without  changing  any  qualitative  results. 

To  reduce  the  order  of  the  system,  we  will  consider  the  limit  of  large  metabolism 
rate,  i.e.  /i  ^  (3,6 ,7.  To  see  how  each  variable  should  scale  with  /1,  we  consider  the 
three-species  equilibrium, 


(H0,So,Po)=(141£s,‘¥1£s). 

At  equilibrium,  there  is  O(g)  more  prey  than  predators.  This  motivates  us  to  scale  the 
prey  variable  in  /i,  so  we  nondimensionalize  the  system  by 

H  1— >  | h  S  1— >  6-s  Pm  ^ X  t  hm  jt. 

The  nondimensional  3D  system  is 

h 

s 

X 

where 

b=  f  >1 

9=~s 

-I- 

We  shall  consider  the  singular  limit  where  a  predator  eats  many  times  in  its  life,  i.e. 
f<Cl.  The  parameter  b  is  the  ratio  of  predator  birthrate  to  death  rate,  g  is  the  ratio 
of  prey  birthrate  to  predator  death  rate,  and  N  is  the  dimensionless  carrying  capacity 
of  the  environment  for  prey.  Note  that  predator  and  prey  populations  are  in  different 
units,  so  their  numerical  values  can  not  be  meaningfully  compared. 

7.2  Reduced  systems 

The  3D  system’s  behavior  can  only  be  well  approximated  by  a  two  dimensional  system 
if  its  behavior  is  roughly  two  dimensional,  for  instance  this  implies  that  it  must  not  be 
chaotic.  There  must  be  a  two-dimensional  slow  manifold  in  phase  space  on  which  all 
solutions  approximately  lie,  possibly  after  some  transient  behavior  as  the  component 
of  the  solution  on  the  fast  manifold  rapidly  decays.  Reduction  of  order  is  achieved  by 


=  —-Xh  +-s  —  h+bs 
e  e 

=  - Xh s  —  s 

e  e 

=  -Xh  +  gX(  l-£), 
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projecting  the  full  system  onto  the  slow  manifold.  To  find  the  slow  manifold,  we  observe 
that  the  system  is  linear  in  h  and  s  if  X  is  regarded  as  a  known  function  of  time: 


The  eigenvalues  of  this  system  are 


A± 


11+1 

e  2 


—  1±1X+1 
e  2 


4  bX 

W+W 


To  first  order  in  e, 

^  =  -1  +  3TTT  +  0<') 


The  A+  eigenvalue  is  0(1),  while  the  A_  eigenvalue  is  negative  and  0(1),  so  solutions 
decay  quickly  along  the  direction  of  the  A_  eigenvector  and  move  more  slowly  in  the 
A+  direction.  Thus  at  a  given  X ,  the  A_  vector  is  tangent  to  the  fast  manifold,  and  the 
A+  vector  is  tangent  to  the  slow  manifold. 


7.2.1  Full  2D  system 

The  A+  eigenvector  yields  a  proportionality  between  s  and  h  on  the  slow  manifold  as  a 
function  of  X, 

S=  2(1  + eft)  (^-1+\/(^  +  l)2+4 ebx)h. 

This  relation  may  be  used  to  reduce  the  dynamical  system  by  eliminating  either  h  or 
s,  but  since  we  are  ultimately  concerned  with  the  total  number  of  predators,  we  define 
Y  =  h  +  s  and  work  in  this  variable.  Applying  the  above  s(h)  relation  to  the  3D  system, 
we  obtain  our  2D  reduction, 


X  =  gX{  1 
Y  =  -Y  + 


_ 2(1  +  eb)XY _ 

N  X  +  1  +  y/(X  +  l)2  +  4  ebX  +  2  eb 
(X  -  1  +  y/(X  +  l)2  +  4ebX)bY 
X  +  1  +  y/(X  +  1)2  +  AebX  +  2eb' 


7.2.2  O(e)  2D  system 


If  we  approximate  the  s(h)  relation  to  O(e),  we  obtain  a  simpler  relation  that  still 
captures  some  affects  of  finite  e, 

s  =  (X  -  ej^x)h  +  0(e2). 

The  O(e)  truncated  reduced  system  is 


X  =  gX(l-  §)- 
bXY 


XY  __  bX3Y 
N>  1  +  X  ~e(X  +  l)3 


Y  =  —Y 


b2X2Y 


1  +  X  (X  +  1)3' 
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7.2.3  Holling  type  II  2D  system 

If  we  truncate  further  and  retain  only  the  0(1)  term  in  e, 

s  =  hX  +  O(e). 


Using  this  simple  proportionality  we  recover  the  Holling  type  II  functional  form. 


1  =  gX(  1  -  § )  - 


Y  =  -Y  + 


bXY 
1  +  X' 


XY 

1  +  X 


The  form  goes  by  various  other  names  in  other  fields,  such  as  the  Jacob-Monod 
form  in  microbiology  or  the  Michaelis-Menten  form  in  enzyme  kinetics. 


7.3  Comparison  of  system  behaviors 

To  understand  what  has  been  lost  by  projecting  onto  the  slow  manifold,  as  well  as 
by  truncating  in  e,  one  must  compare  the  behavior  of  the  full  3D  system  with  the 
behaviors  of  the  reduced  2D  systems.  There  is  no  unique  measure  of  the  quality  of  the 
approximation;  different  properties  are  approximated  better  than  others,  so  the  value 
of  the  approximation  depends  ultimately  on  what  properties  are  of  interest. 

For  the  3D  system,  we  ignore  the  separate  dynamics  of  the  h  and  s  variables  and 
consider  only  their  sum,  Y,  because  this  is  the  quantity  of  interest  and  the  one  that 
compares  directly  to  the  2D  models.  As  we  demonstrate  below,  all  four  models  have 
the  same  qualitative  behavior,  some  representative  phase  portraits  of  which  are  given  in 
Figure  7.  There  is  a  trivial  equilibrium  at  the  origin  representing  mass  extinction,  but 
it  is  always  unstable.  There  is  a  prey-only  equilibrium,  stable  only  when  b  and  N  are 
small.  That  is,  when  predators  are  not  born  too  fast,  and  saturation  population  of  prey 
is  not  too  large,  both  factors  that  would  inhibit  predator  success.  When  b  and  N  are  a 
bit  larger  so  conditions  are  a  bit  better  for  predators,  the  prey-only  equilibrium  becomes 
unstable  as  a  stable  two-species  equilibirum  becomes  physical  and  splits  off  from  it  in 
a  transcritical  bifurcation.  The  two-species  equilibrium  is  a  plain  sink  initially  and 
becomes  a  spiral  sink  at  larger  b  and  N.  When  b  and  N  are  increased  further  still, 
the  two-species  equilibrium  undergoes  a  Hopf  bifurcation,  losing  its  stability  to  a  limit 
cycle.  All  orbits  are  bounded  for  all  parameter  values.  Since  all  four  systems  share 
this  qualitative  picture,  the  effects  of  approximation  appear  only  in  the  quantitative 
differences  between,  say,  equilibria  and  limit  cycle  locations,  or  bifurcation  values. 


7.3.1  Equilibria 

At  all  parameter  values,  all  four  systems  have  equilibria  at  the  origin  and  at  (X,  Y)  = 
( N ,  0).  Solving  for  the  nontrivial  two-species  equilibria,  we  obtain  Yq  =  r^j(l  —  for 
all  four  systems,  though  the  Xo  value  may  vary  between  systems.  Clearly  the  nontrivial 
equilibrium  is  only  physical  (Ij)  >  0)  when  X0  is  less  than  N,  the  prey-only  saturation 
population.  That  is,  such  equilibria  never  represent  mutualistic  solutions.  The  Xo 
values  for  the  two-species  equilibria  in  each  system  are  tabulated  in  Table  2.  The  Xo 
value  for  the  0(e)  system  is  more  cleanly  expressed  implicitly  by  the  cubic  equation 
of  which  it  is  the  only  positive  real  root.  The  full  2D  system  has  exactly  the  same 
equilibrium  as  the  3D  system,  while  the  truncated  systems  have  different  two-species 
equilibria,  which  converge  to  the  3D  value  as  e  — >  0. 
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Figure  7:  Sample  phase  portraits  representing  all  four  possible  qualitative  behaviors,  con¬ 
sidering  N  as  the  control  parameter:  stable  prey-only  equilibrium  (upper  left),  stable  two- 
species  equilibrium  with  plain  sink  (upper  right),  stable  two-species  equilibrium  with  spiral 
sink  (lower  left),  and  a  limit  cycle  around  the  two-species  equilibrium  (lower  right).  These 
portraits  were  generated  by  integrating  the  3D  system  with  b  =  2,  g  =  1,  and  e  =  0.5. 
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Table  2:  Locations  of  the  two-species  equilibria  in  the  four  models. 


Model 

*0 

Y0 

3D 

1+e 

6—1 

bg  ( 1  l+i 

6— 1  v  N(b—1)  ) 

Full  2D 

l+e 

6-1 

bg  ( 1  1+7  L- 

6-1  l1  TV  (6 — 1)3 

0(e)  2D 

(b  -  1)X^  +  (26  -  eb2  -  3)Xq  +  (6  -  3)X0  -1  =  0 

bg  (1  XD 

6-1  N  ) 

Holling  type  II  2D 

1 

6-1 

bg  (1  i  T” 

6—1  V1  TV  (6 — 1)3 

7.3.2  Linear  stability 

To  analyze  the  linear  stability  of  the  three  equilibria  in  all  four  systems,  we  linearize 
each  system  about  an  arbitrary  point,  (ho,  Sq,  X0),  or  (X0,lo). 


3D  system 


i  (S) ' 
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About  the  origin,  every  linearized  system  has  an  eigenvalue  of  g ,  which  is  positive, 
so  the  extinction  equilibrium  is  always  unstable.  Linearizing  about  the  prey-only  equi¬ 
librium,  one  finds  that  each  system  goes  unstable  according  to  the  same  parameter 
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inequalities  that  determine  when  the  two-species  equilibrium  exists.  This  is  as  expected 
from  a  transcritical  bifurcation;  the  stability  of  the  prey-only  equilibrium  changes  pre¬ 
cisely  when  it  collides  with  the  two-species  equilibrium,  which  is  also  the  moment  when 
that  equilibrium  becomes  physical.  As  for  the  two-species  equilibrium,  we  can  study  its 
stability  analytically  in  the  Holling  type  II  system,  but  the  other  linearized  systems  are 
messy,  so  we  solve  their  stability  eigenproblems  numerically. 


7.3.3  Stability  of  the  two-species  equilibrium  in  the  Holling  type  II 
2D  system 


The  Holling  type  II  system  linearized  about  its  two-species  equilibrium  is 


d_ 

dt 


(“) 


'  3  (\ _ b+1 _ 3 

H  N(b-l)J 


The  stability  of  the  two-species  equilibria  is  not  hard  to  compute  analytically  for  the 
Holling  type  II  system.  The  characteristic  equation  of  the  linearized  system  is 

A2  “  b  (1  ~  Nb(b-1)  )  A+f(5_1_F)=°- 

For  the  two-species  equilibrium  to  exist,  the  0(1)  coefficient  of  the  characterisitc  equa¬ 
tion  must  be  positive.  Thus,  solving  the  quadratic  equation  for  A,  the  discriminant  will 
either  be  imaginary  or  of  smaller  magnitude  than  the  O(A)  coefficient.  Either  way,  both 
eigenvalues  will  be  negative  (the  equilibrium  will  be  stable)  if  and  only  if  the  0(A)  co¬ 
efficient  is  positive,  i.e.  when  N  <  When  N  exceeds  this  value,  all  three  equilibria 
are  unstable.  We  later  prove  that  all  orbits  are  bounded,  so  the  Poincare-Bendixon 
theorem  will  guarantee  that  the  system  converges  to  a  limit  cycle. 


7.3.4  Bifurcations 

The  transcritical  bifurcation  occurs  when  Yq  exceeds  zero,  which  occurs  when  the  Xo 
expressions  reported  in  Table  2  are  less  than  N.  For  each  system,  this  is  possible  only 
when  b  >  1.  The  exact  relation  between  b  and  N  at  the  bifurcations  are  given  in  Table 
3.  Note  that  from  a  point  in  parameter  space  where  the  prey-only  equilibrium  is  stable, 
the  bifurcation  may  be  produced  by  increasing  either  b  or  N.  The  point  in  parameter 
space  where  the  Hopf  bifurcation  occurs  has  a  simple  analytic  expression  for  the  Holling 
type  II  model,  so  this  is  also  given  in  Table  3.  For  the  other  models,  the  N  at  which  the 
Hopf  bifurcation  occurs  was  computed  numerically  for  given  values  of  6,  g  and  e,  and 
some  representative  results  are  plotted  in  Figure  8.  The  3D  transcritical  bifurcation 
depends  only  on  6,  N  and  e,  while  the  Hopf  bifurcation  depends  also  on  g,  but  quite 
weakly  so.  It  is  clear  from  Figure  8  that  the  full  2D  model  captures  the  transcritical 
bifurcation  perfectly,  while  the  truncated  models  are  inaccurate  when  e  becomes  large. 
At  the  Hopf  bifurcation,  the  full  2D  model  captures  the  3D  behavior  imperfectly,  but 
again  much  better  than  the  truncated  models. 


7.3.5  Lyapunov  stability 

Each  system  undergoes  only  the  two  bifurcations  we  have  studied  and  has  no  other 
fixed  points.  All  that  remains  is  to  verify  that  in  each  system  the  orbits  are  bounded 
for  all  parameter  values.  We  shall  do  this  by  the  Lyapunov  method  for  the  3D  system 
and  the  Holling  type  II  2D  system.  We  shall  not  prove  boundedness  for  the  other  two 
systems,  whose  algebraic  nonlinearities  would  make  it  a  cumbersome  task,  but  we  can 
feel  confident  in  its  veracity. 
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Table  3:  Bifurcation  points  of  the  four  models.  For  the  0(e)  model,  Xq  is  defined  implicitly 
by  the  formula  given  in  Table  2. 


Model 

Transcritical  bifurcation 

Hopf  bifurcation 

3D 

0  <  jgf  <  N 

Found  numerically 

Full  2D 

0  <  <N 

O(e)  2D 

0  <  Xo  <  N 

Holling  type  II  2D 

«<X<N 

gfXv 

Holling  type  II  2D  system  Examining  the  X  equation,  we  see  that  X  <  0 
whenever  X  >  N.  Thus,  if  X  <  IV  at  the  initial  condition,  it  remains  true  for  all  time. 
To  put  an  upper  bound  on  Y  that  is  valid  for  all  parameters,  we  must  consider  X  and 
Y  together  in  a  Lyapunov  functional. 

Let  L  =  X  +  jY .  The  proportionality  between  X  and  Y  is  chosen  such  that  the 
nonlinear  Holling  type  II  terms  in  L  cancel: 

L  =  gX(l-§)-lY 


Our  goal  is  to  bound  L  by  an  affine  function  of  L,  i.e.  L  <  a  —  (3L,  where  (3  >  0.  This 
will  imply  that  L  <  %  for  all  time  if  it  is  true  initially.  To  bound  L  by  such  a  term  we 
must  bound  the  quadratic  X  term  by  an  affine  function  of  X  with  a  negative  coefficient 
on  X.  Any  line  tangent  to  the  parabola  at  X*  >  ^  will  suffice,  but  we  seek  the  line 
that  minimizes  thereby  providing  the  optimal  bound  on  L.  An  arbitrary  tangent 
line  gives  the  bound 


which  induces  a  bound  on  L, 


L  ^  Tv 


-M}*- 


Thus, 


r  •  r  X* 

L  <  mt  - - ,  ,  . 

x,  >n/2  N  min  {  -  1,1} 


Assuming  the  above  infonum  occurs  at  an  X*  such  that  —  1  <  \ ,  the  optimal  choice 
of  X*  in  fact  contradicts  the  assumption  when  b  >  1.  So,  the  minimum  must  be 
meaning  that  X*  >  X  (f  _|_  I)  >  jy.  The  bound  on  L  then  becomes 


L  <  inf 

x,>f(i+L 


=  K(b. 

N  4  b\u 


l)2 


Putting  this  in  terms  of  X  and  Y,  and  adding  the  known  bound  on  X  alone, 

X  <  min  j  jpdl  +  b)2  —  \Y,N |. 


3D  system  Let  L  =  h  +  as  +  f3X.  To  make  the  hX  terms  to  vanish  from  L,  we  let 
/3  =  which  clearly  requires  a  >  1  for  (3  to  be  positive.  Thus, 


L  =  —  [T  +  1  —  (  \  +  b) / a\as  —  h+  -(a  —  l)gX(l  -  y) 
<  -m(a)L+  i(a-  l)X[(g  +  m{aj)  -  ^] , 
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Figure  8:  Values  of  N  as  a  function  of  e  at  the  transcritical  and  Hopf  bifurcations  for  all 
four  systems  with  6  =  2  and  g  =  1.  For  the  transcritical  bifurcation,  the  full  2D  systems 
coincides  with  the  3D  system. 


where  m(a)  =  min  {1,^  +  1  —  (-  +  6) /a}.  We  now  bound  the  quadratic  X  term  by  its 
maximum  value  without  bothering  to  find  the  optimal  bound  linear  in  X. 

L  <  -m{a)L+  1  )(g  +  m{a))2. 

This  yields  a  bound  on  L  that  can  be  optimized  over  all  allowable  a: 

L<  inf  ~r~  a7\  (g  +  m(a))~ . 
a>l  4£9  m(“)  v  '> 


Assuming  the  infimum  occurs  when  m(a)  =  1  implies  that  the  optimal  bound  is  obtained 
by  choosing  a  =  1+,  which  contradicts  m(a )  =  1.  Let  a  >  giving  our  best  result 

for  a  bound  on  the  Lyapunov  functional: 


L  < 


inf 


a> 


l  +  eb 

l  +  « 


N  a- 1 
4  e2g  1  +  e  - 

Oi 


(g+l  +  e 


l+eb')2 

Ct  ) 


The  optimal  bound  can  be  calculated  given  the  other  parameters,  but  the  more  impor¬ 
tant  conclusion  is  that  some  such  finite  bound  always  exists  for  L. 


7.3.6  Summary  of  results 

Our  analysis  strongly  suggests  that  X  and  Y  have  the  same  qualitative  behavior  in  all 
four  systems,  though  to  make  this  result  rigorous,  one  needs  Lyapunov  bounds  on  the 
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full  2D  and  O(e)  2D  systems,  and  one  needs  to  show  analytically  that  a  Hopf  bifurcation 
occurs  in  all  systems  as  it  does  in  the  Holling  type  II  system.  Trusting  that  the  systems 
indeed  all  have  the  same  behavior,  they  differ  only  quantitatively.  We  have  seen  that 
the  full  2D  system  has  the  same  two-species  equilibrium  as  the  3D  system,  while  the 
truncated  systems  do  not,  and  quantitative  differences  in  bifurcation  values  were  shown 
already  in  Figure  8. 

Phase  portraits  produced  by  the  different  models  appear  in  Figure  9.  Although  the 
2D  systems  began  with  the  same  initial  conditions,  we  must  compare  them  each  to 
their  own  corresponding  3D  solution  because  they  each  correspond  to  slightly  different 
decompositions  of  Y0  into  ho  and  So-  However,  the  three  3D  solutions  tend  to  be  quite 
similar.  The  top  row  of  Figure  9  shows  solutions  at  small  N,  when  the  prey-only 
equilibrium  is  stable.  The  prey-only  equilibrium  is  identical  in  all  four  systems,  so  the 
phase  portraits  agree  well  even  for  e  =  0.5.  The  middle  row  of  Figure  9  shows  solutions 
for  larger  N,  when  the  two-species  equilibrium  is  stable.  The  different  systems  agree 
well  when  e  is  0.05,  but  at  0.5  the  locations  of  the  equilibria  differ  significantly,  so  the 
respective  phase  trajectories  spiraling  towards  them  are  quantitatively  quite  different. 
The  bottom  row  of  Figure  9  shows  long-time  solutions  at  still  larger  N.  When  e  is  0.05, 
the  limit  cycles  of  the  full  2D  and  0(e)  2D  systems  approximate  the  3D  limit  cycle 
quite  well,  while  the  Holling  type  II  system  does  a  bit  worse.  When  e  is  0.5,  the  Holling 
type  II  system’s  limit  cycle  is  much  too  large,  the  full  2D  system’s  is  too  small  but  a 
bit  better,  and  the  O(e)  2D  system  has  not  yet  gone  through  the  Hopf  bifurcation. 

The  full  2D  system  approximates  the  3D  system  well  for  e  <  0.5,  while  the  Holling 
type  II  2D  system  is  quantitatively  accurate  only  when  e  is  an  order  of  magnitude 
smaller.  We  wish  to  extrapolate  these  truths  to  other  models  where  the  Holling  type 
II  or  full-order-in-e  functional  forms  might  be  used  as  predation  laws  without  repeating 
their  rigorous  derivation  from  a  higher-order  dynamical  system.  If  e  is  very  small, 
or  if  one  is  only  concerned  with  qualitative  features,  as  is  often  the  case  in  biological 
modeling,  the  Holling  type  II  functional  form  is  certainly  satisfactory.  If  e  is  closer 
to  unity,  and  the  quantitative  properties  of  the  system  matter,  as  is  often  the  case  in 
enzyme  kinetics,  the  full-order-in-e  functional  form  would  be  a  better  choice.  The  O(e) 
functional  form  probably  offers  neither  enough  simplicity  nor  accuracy  to  be  chosen 
over  the  other  two. 
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e  =  0.05,  N  =  5 


X  (prey) 


Figure  9:  Phase  portraits  of  all  four  systems  for  values  of  N  with  e  =  0.05  (left)  and  e  =  0.5 
(right),  and  b  =  2,  g  =  1.  All  solutions  began  at  (2,2),  but  only  the  late-time  behavior 
is  shown  in  the  bottom  two  plots  to  make  the  limit  cycles  clear.  The  asterisks  are  the 
equilibria  of  the  different  systems,  which  always  coincide  for  the  3D  system  and  the  full  2D 
system. 
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